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ABSTRACT 

The application of orthonormal basis functions such as 
Prolate Spheroidal Wave Functions (PSWF) for accurate 
source modeling in radio astronomy has been compre- 
hensively studied. They are of great importance for high 
fidelity, high dynamic range imaging with new radio tele- 
scopes as well as conventional ones. But the construction 
of PSWF is computationally expensive compared to other 
closed form basis functions. In this paper, we suggest a 
solution to reduce its computational cost by more eflftcient 
construction of the matrix kernel which relates the image 
domain to visibility (or Fourier) domain. Radio astronom- 
ical images are mostly represented using a regular grid of 
rectangular pixels. This is required for efficient storage and 
display purposes and moreover, comes naturally as a by 
product of the Fast Fourier Transform (FFT) in imaging. 
We propose the use of Delaunay triangulation as opposed 
to regular gridding of an image for a finer selection of the 
region of interest (signal support) during the PSWF kernel 
construction. We show that the computational efficiency 
improves without loss of information. Once the PSWF 
basis is constructed using the irregular grid, we revert 
back to the regular grid by interpolation and thereafter, 
conventional imaging techniques can be applied. 

Index Terms- Radio astronomy. Radio interferometry, 
Deconvolution 

I. INTRODUCTION 

For radio interferometric imaging, we measure the spa- 
tial coherency or Fourier components of the celestial 
sources on diff'erent baselines in our array (also called 
as visibilities). Due to a finite number of receivers, we 
only sample the Fourier components at a finite number of 
positions. After calibration or correction for the corruptions 
due to the propagation path and instrumental eff'ects, we 
have to deconvolve the source flux from the antennae 
response pattern to obtain its true intensity map or its radio 
image. 
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Conventionally, CLEAN fll,(2l algorithm has been used 
for modeling and removal of the known sources such that 
the faint unknown sources in the background can appear. 
It is a pixel based deconvolution approach. This means 
any point in the image is associated with a weighted 
delta function. CLEAN will fail to perform satisfactorily, 
if two point sources are too close or if the source is 
partially resolved, as it has analytically been proven in O. 
Moreover, we have shown that making image grid smaller 
(expecting that accuracy of modeling extended sources 
improves), is a futile eflfort by presenting the Cramer Rao 
Bound (CRB) analysis |3|. Instead, we showed that using 
a 2-D orthonormal basis for deconvolution achieves the 
theoretical noise bound or CRB, regardless of the source 
structure 0. 

Cartesian shapelets have been used for source modeling 
in 1 3 1, (and references therein). It has been proven that 
shapelets provide a way to model extended sources by a 
closed form construction of an orthonormal basis. But this 
technique requires additional constraints such as informa- 
tion theoretic bounds to limit the model order. In IH, 0, 
we have shown that we can construct the optimal basis 
to model the source by minimum number of basis func- 
tions and less artifacts outside Region Of Interest (ROI) 
using prolate spheroidal wave functions. In addition, unlike 
shapelets, they do not require additional explicit constrains 
from information theory, as it is implicitly considered 
in their mathematical derivation (see |5|). However, the 
construction of PSWF is computationally more expensive 
as compared to shapelets (H. 

In this paper, we present computationally efl&cient con- 
struction of PSWF suitable for modeling extended sources. 
We achieve this by downsampling both the Fourier plane 
visibilities as well as image pixels. Most astronomical (as 
well as other) images are represented by a regular 2-D 
array of rectangular pixels. In contrast, we propose to use 
Delaunay triangulation to construct an irregular grid with 



fewer number of pixels outside ROI and finer selection 
of the image inside ROI to represent the image during 
the construction of PSWF. We can control the size of the 
triangulation without losing information. Once the PSWF 
basis has been constructed in this manner, we revert back 
to the original pixel grid by interpolation. 

Notation: We denote vectors in bold lowercase and ma- 
trices in bold uppercase. The matrix transpose, Hermitian, 
and pseudoinverse are denoted by (.)^ , (.)^ and (.)^ 
respectively. Operation ||.|| gives the Frobenious norm of 
a matrix. 

II. MATHEMATICAL PRELIMINARIES 
II-A. Radio interferometric imaging 





(a) 



(b) 



Fig. 1. (a) Sampling points (total Na) in the Fourier 
(visibility) plane, (b) Image of N rectangular pixels, with 
the support (ROI) area in white. The ROI area has Nt 
pixels. 

In radio interferometric imaging (see Fig.[T]), the relation 
between the intensity of the ^-th pixel, /(/^,m^) in an 
image and the p-th point, f(up, Vp) in the visibility domain 
is described by van Cittert-Zernike theorem [6J as: 

Na-l 

f(l,, m,) = £ f(up, v^)^^-2-(V.--.v,)^ 

N-l 

f(up,vp) = ^/(/„m,)^-^'2^(^^"^"-^^^>, 

pe[0,Na-ll qe[0,Nl 

where A^^^ is the total sampling points in visibility domain. 
N is the total number of pixels in the image. If we rewrite 
them in vectorized form, we will obtain: 



(2) 



f - Tf , f 

f = [/(/o,mo), ...,/(/A^-i,mA^_i)]^ 

The matrix T (size Na x N) contains ^-727r(/^W/^+'^?V;,) 
its ^-th row and p-th column. 



II-B. Prolate spheroidal basis 

For a complete mathematical derivation of the PSWF, 
the reader is referred to |31 and references therein. In 
the following, we will mention a few key mathematical 
relations of PSWF derivation for stating our problem. 

We aim to obtain representative basis, p(l,m) for a 
source that maximizes the energy in our ROI i.e. within the 
boundary by which a source is recognized in an image. The 
larger the ratio given in ([3]), the less artifacts will remain 
outside the ROI: 



m)eROI 



\p(l, m) 



The vectorized form of which can be written as: 



(3) 



(4) 



where ij, (size x A/^) is a selection matrix to select the 
pixels that belong to the ROI. is the number of pixels in 
the ROI. Then the PSWF basis can be obtained by eigen- 
decomposition of the kernel, K (size Nb x Nt) given as: 



K = irT^(TT^)"^TI/, 



(5) 



Most of the computations are required first to find the 
pseudo-inverse term, (TT^)^ (size A^^ x A^^) and secondly, 
to find the eigenvalue decomposition of K (size N^xNt). 
Note that the rank of (TT^)^ is at most N and usually A^^ » 
N > Nt.By properly downsampling both the Fourier plane 
sampling points (reducing Na) as well as the image pixels 
(reducing N and A^^) we can significantly cut down the 
cost of computations. In the next section, we will explain 
how it has efficiently been done, followed by an example. 

III. COMPUTATIONAL COST REDUCTION 

We take several steps to reduce the computational cost 
of ([5]). The main criterion for this reduction is given by the 
Landau-Pollak theorem |7|. In its simplest form (applicable 
to radio interferometry) [4J, Landau-Pollak theorem states 
that the number of degrees of freedom of any given source 
with compact support. No could be no greater than the 
product of area of the source. Aim and the area of the source 
observed in the Fourier plane, A^v i.e. No < Aim x A^v- 

III-A. Reducing A^^ 

Since we already have an irregular set of A^^ sampling 
points on the Fourier plane, it is straightforward to down- 
sample this to a lower value A^^, which is already presented 
in |5|. However, we can do even better by combining this 
with imaging weights. Each point on the Fourier plane, 
(say p) has a weight associated with it, p(up,Vp) e [0, 1]. 




Prior to taking the Fourier transform in ([T]), we multi- 
ply each data point f(up,Vp) by this weight. When we 
downsample, we generate a uniformly distributed random 
number r(up,Vp) e [0,1] and if r(up,Vp) < p(up,Vp) we 
select the p-th point. 

III-B. Reducing N but finer selection of ROI 

The motivation behind reducing the number of image 
pixels is that not all the pixels are required to define 
the ROI (because of Landau-Pollak theorem) while we 
construct the PSWF. In Fig. [2] (a), we have shown the ROI 
of the source as well as the Point Spread Function (PSF). 
The amount of information can be loosely translated to 
the number of PSFs that we can pack within the ROI, as 
shown in Fig. [2] (b). The optimal way to pack the PSFs 
within the ROI boils down to a circle packing problem, 
which is generally understood as NP-hard. 
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Fig. 2. (a) The ROI of an extended source and the PSF 
(the ellipse), (b) The ROI covered by the PSF (an ellipse 
packing). 

Therefore, we follow a heuristic approach: Instead of 
solving a packing problem, we construct an irregular grid 
of triangles to cover the ROI as well as the area of the 
image outside the ROI. Obviously, we employ Delaunay 
triangulation |8| for this purpose. We shall illustrate this 
with an example as follows: In Fig. [51 we have shown an 
image of an extended radio source. This image (dimensions 
118x64) has 7552 pixels. 

In Fig.m we have shown the ROI for this source, which 
has 1826 pixels to represent the amount of information 
present. Therefore, we have generated the Delaunay trian- 
gulation for the ROI and its exterior as shown in Fig. \5\ 
such that the triangle size within the ROI is smaller than 
the exterior triangles. We take the centroids of the triangles 
as our reduced grid to generate the PSWF. Therefore, we 
end up with a total number of pixels of 1220 while the 
ROI has 1080 pixels. 

III-C. Reducing cost of pseudoinverse 

The cost of dll) mostly involves evaluation of the pseu- 
doinverse (TT^)^ Instead of directly evaluating this, we 
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Fig. 3. A test image of an extended radio source, with 
dimensions 118 x 64 pixels. 
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Fig. 4. The ROI of the source in Fig. [3] The total image 
has 7552 pixels while the ROI has 1826 pixels. 
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Fig. 5. Delaunay triangulation of the image in Fig. [3l 
The triangulation inside the ROI has 1080 triangles while 
the triangulation outside the ROI has 140 triangles. The 
triangle size within the ROI is made smaller than the 
exterior triangles. 



solve a set of linear equations such as 

(TT^ + rl)x = Ax = y (6) 

where y is a small positive regularization parameter. De- 
pending on the context, we are given y and we need to find 
X. The solution of the linear system for x is done using 
Q-less QR factorization. 

Given QR = A (where Q is unitary and R is upper 
triangular), we first solve 

R^Rx = A^y (7) 

for X and make one step correction 



(or values of a and b). If the basis vectors are orthonormal, 
we will obtain the lowest value for the Frobenious norm. 
Fig. 7 (a) shows that the optimality of the source modeling 
is not aff'ected significantly from one setting to another, 
because the basis vectors are orthonormal with negligible 
errors. Thus, the CRB remains minimum |3 1. Moreover, the 
condition number of P^P helps us to understand how its 
condition aff'ects the further computation. Ultimately, these 
conclude that the setting with a = 0.2 and /? = 2.1 is the 
most suitable one to fulfill the aforementioned criteria. By 
this setting, we downsample to N' = 3270 pixels in which 
the ROI has A^^ = 3199 pixels. In contrast, the original 
pixel grid has = 7552 and Nt = 1826. 



r = y - Ax 
solve R^R'e = A^r 



(8) 



X = X + e 



to get the solution x. 



IV. RESULTS 

To demonstrate how our solution works in practice, we 
give results on the construction of PSWF for the extended 
source shown in Fig. 3. The source was observed by 
the LOFAR (http://www.lofar.org) radio telescope during 
system testing. For an 8 hours of observation, we obtain, 
Na = 7896636 data points on the Fourier plane. We reduce 
this to N'^ = 11041 points by applying our downsampling 
scheme as explained in section 3.1. 

We use |8| to generate conforming Delaunay triangula- 
tions for the ROI and the exterior of the image. Selection 
of triangle scales for inside and outside of the ROI (see 
Fig.5), has to be done in such a way that the least amount 
of information is lost. For this observation, the PSF minor 
diameter is about d = 4.2x 10"^ radians. To achieve the 
most appropriate triangulation, we have set diff'erent scales 
for the triangles inside the ROI (a x d) and outside the ROI 
(bxd), where a,b are the parameters that we vary. 

We show the eigenvalue spectrum of the kernel K, for 
diff'erent values of a and b in Fig. 6. The eigenvalue 
spectrum gives us a measure of information a certain basis 
can represent. By the results shown in Fig. 6, we see that in 
most settings, we get about 100 basis functions that carry 
significant information, which is close to the limit given 
by Landau-Pollak theorem for this case. 

Once we obtain the basis functions for the downsampled 
grid, we interpolate back to the original image grid to 
get the basis vectors corresponding to the original image. 
Let P be the matrix whose columns correspond to the 
interpolated basis vectors. Then, we study the Frobenious 
norm of P^P-I and the condition number of P^P in Fig. 7. 
This, in addition to the eigenvalue spectrum in Fig. 6 give 
us a benchmark to choose the most suitable triangulation 
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Fig. 6. Eigenvalue spectrums for diff'erent triangulations, 
with varying a and b. 



V. CONCLUSIONS 

The use of PSWF for source modeling and deconvolu- 
tion of extended sources was proposed earlier in 0| and 
|5|. In this paper, we suggested a scheme to reduce its 
computational cost. By using the Delaunay triangulation, 
we can significantly reduce the cost of construction of 
PSWF. We have also shown by a real example that we 
do not lose information by this cost reduction. Further 
improvement will focus on finer selection of the ROI 
depending on the source structure such that an area with 
higher intensity is represented with more triangle pixels. 

VI. REFERENCES 

[1] J. A. Hogbom, "Aperture synthesis with a non-regular 
distribution of interferometer baselines," A6rA SuppL, 
vol. 15, pp. 417-426, Jun. 1974. 





Fig. 7. (a) Frobenious norm of P^P - I (scaled by 1 /N^^^.^ where Ntasis is the number of basis functions) for different 
triangulations (varying a and b) of the image, (b) Condition number of P^P for different triangulations. 



[2] U. Schwarz, "mathematical statistical description 
of the iterative beam removing technique (method 
clean)," A&A, vol. 65, no.l, pp. 345-356, 1978. 

[3] S. Yatawatta, "Fundamental limitations of pixel based 
image deconvolution in radio astronomy," in proc. 
IEEE Sensor Array and Multichannel Signal Process- 
ing Workshop (SAM), Jerusalem, Israel, pp. 69-72, 
Nov. 2010. 

[4] , "Shapelets and related techniques in radio astro- 
nomical imaging," in proc. URSI GA 2011, Istanbul, 
Turkey, Aug. 2011. 

[5] , "Radio astronomical image deconvolution using 

prolate spheroidal wave functions," in proc. IEEE 
International Conference on Image Processing (ICIP), 
Brussels, Belgium, pp. 2813-2816, Sep. 2011. 

[6] W. N. Brouw, "Aperture synthesize," Methods in Com- 
putational Physics, vol. 14, pp. 131-175, 1975. 

[7] H. J. Landau and H. O. Pollak, "Prolate spheroidal 
wave functions, fourier analysis, and uncertainty-ii," 
Bell Syst. Tech. J., vol. 40, no. 2, p. 6584, 1961. 

[8] P. Persson and G. Strang, "A simple mesh generator 
in matlab," SIAM Review, vol. 46, no. 2, 2004. 



